#include "device.h"

/* -------------------------------------------- */
/* -------------------------------------------- */

void Device::q_eraser() { /* Destroy the quantum structure */
  int j,i,k;
  if(ynum_mesh_points<=1 && tnum_mesh_points==0) {
    for(k=0;k<num_lyr;k++) {
      for(i=0;i<lyr[k].numx();i++) {
	if(lyr[k].grid1d[i].qptr_e!=NULL) {
	  delete lyr[k].grid1d[i].qptr_e;
	  lyr[k].grid1d[i].qptr_e = 0;
	}
	if(lyr[k].grid1d[i].qptr_h!=NULL) {
	  delete lyr[k].grid1d[i].qptr_h;
	  lyr[k].grid1d[i].qptr_h = 0;
	}
      }
    }
  }
  else if(ynum_mesh_points>1 && tnum_mesh_points==0) {
    for(k=0;k<num_lyr;k++) {
      for(i=0;i<lyr[k].numx();i++) {
	for(j=0;j<ynum_mesh_points;j++) {
	  if(lyr[k].grid2d[i][j].qptr_e!=NULL) {
	    delete lyr[k].grid2d[i][j].qptr_e;
	    lyr[k].grid2d[i][j].qptr_e = 0;
	  }
	  if(lyr[k].grid2d[i][j].qptr_h!=NULL) {
	    delete lyr[k].grid2d[i][j].qptr_h;
	    lyr[k].grid2d[i][j].qptr_h = 0;
	  }
	}
      }
    }
  }
  else nrerror("Q_ERASER: 3D GRID NOT IMPLEMENTED YET!");
    
}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* FIND_HETERO_CB returns "left[n][2]" and "right[n][2]" boundaries of "n" quantum structures
   presented in the conduction band of the device. For quantum structure we mean a single block where 
   Schrodinger equation has to be solved wholly. 
   It first spots where possibly heterointerfaces are. 
   By changing the parameter secure_length, we modify the width of the distance of the left or right
   boundary from the real quantum structure. 
   The "max_length" instead checks whether to include the next interface in the 
   quantum structure. ATTENTION: max_length>secure_length!
   This is done to avoid neglecting tails of wave functions.
   The pointer "fp" links to the debugging purpose file "quantum.log" */

int Device::find_hetero_cb(int ** left, int ** right, FILE * fp) {

  const int MAX_STEPS = 10;
  const double band_step = 0.2;  // eV

  int i,k,j,n,l,p,real_min[2],min[MAX_STEPS][2],lmin[MAX_STEPS][2],rmin[MAX_STEPS][2];
  double d,cb1,cb2,min_val;

  double max_length = sim->Max_quantum_width; /* cm */
  double secure_length = sim->Secure_length; /* cm */


  /* Find all the steps in the band structure */
  n = 0;
  cb1 = lyr[0].grid1d[0].chi();
  for(k=0;k<num_lyr;k++) {
    for(i=0;i<lyr[k].numx();i++) {
      cb2 = lyr[k].grid1d[i].chi();
      if( fabs(cb1-cb2) > band_step ) {
	min[n][0]=k;
	min[n][1]=i;
	n++;
      }
      cb1 = cb2;
    }
    if(n>=MAX_STEPS) {
      printf("FIND_HETERO_CB: Too many heterostructures\n");
      exit(1);
    }
  }

  if(n==0) {
    if(fp) fprintf(fp,"FIND_HETERO_CB: No heterostructures\n");
    return 0;
  }

  if(fp) printf("Num step = %d in %d,%d\n",n,min[n-1][0],min[n-1][1]);

  j = 0;
  /* LOOP ON STEPS */
  while(j<=n-1) {
    /* Set the left boundary position */
    if(min[j][1]==0) {
      lmin[j][0]= min[j][0]-1;
      lmin[j][1]= lyr[min[j][0]-1].numx()-1;
    }
    else {
      lmin[j][0]= min[j][0];
      lmin[j][1]= min[j][1]-1;
    }
    d = lyr[lmin[j][0]].grid1d[lmin[j][1]].xsize();
    while( d<secure_length && 
	   ( (j==0 && (lmin[j][0]>0 || (lmin[j][0]==0  && lmin[j][1]>0))) || 
	     (j>0 && ( lmin[j][0]>rmin[j-1][0] || (lmin[j][0]==rmin[j-1][0] && lmin[j][1]>rmin[j-1][1]) ) ) 
	   ) 
	 ) {
      if(lmin[j][1]==0) {
	lmin[j][0] = lmin[j][0]-1;
	lmin[j][1] = lyr[lmin[j][0]-1].numx()-1;
      }
      else lmin[j][1] = lmin[j][1]-1;
      d += lyr[lmin[j][0]].grid1d[lmin[j][1]].xsize();
    };
    if(lmin[j][0]==0 && lmin[j][1]==0) nrerror("FIND_HETERO_CB: The quantum structure has reached the left edge of the device. Decrease the secure_length variable or put an extra layer at the left of the device.");

    /* Check if within a max_length there is another step and change right boundary */
    rmin[j][0] = min[j][0];
    rmin[j][1] = min[j][1];
    k = min[j][0];
    i = min[j][1];
    d = 0.5*lyr[k].grid1d[i].xsize();
    if(i<lyr[k].numx()-1) i++;
    else {k++; i=0;}
    p = 1;
    while( d<max_length && k<num_lyr) {
      d += lyr[k].grid1d[i].xsize();
      if(i<lyr[k].numx()-1) i++;
      else {k++; i=0;}
      if(i==min[j+p][1] && k==min[j+p][0]) {
	rmin[j][0] = min[j+p][0];
	rmin[j][1] = min[j+p][1];
	d = 0;
	for(l=j+p;l<n-1;l++) {
	  min[l][0]=min[l+1][0];
	  min[l][1]=min[l+1][1];
	}
	n--;
      }
    }
  
    /* Set the right boundary position */
    d = 0;
    d += lyr[rmin[j][0]].grid1d[rmin[j][1]].xsize();

    while(d<secure_length && (rmin[j][0]<num_lyr-1 || (rmin[j][0]==num_lyr-1 && rmin[j][1]<lyr[num_lyr-1].numx()-1))) {
      if(rmin[j][1]==lyr[rmin[j][0]].numx()-1) { 
	rmin[j][0] ++;
	rmin[j][1] = 0;
      }
      else rmin[j][1]++;
      d += lyr[rmin[j][0]].grid1d[rmin[j][1]].xsize();
    }
    if(rmin[j][0]==num_lyr-1 && rmin[j][1]==lyr[num_lyr-1].numx()-1)   
	nrerror("FIND_HETERO_CB: The quantum structure has reached the right edge of the device.Decrease the secure_length variable or put an extra layer at the right of the device.");

    /* Check whether there is a minimum */
    k = lmin[j][0];
    i = lmin[j][1];
    real_min[0]=k;
    real_min[1]=i;
    min_val = lyr[k].grid1d[i].cb();
    while(k<rmin[j][0] || (k==rmin[j][0] && i<=rmin[j][1])) {
      if(lyr[k].grid1d[i].cb()<min_val) {
	real_min[0]=k;
	real_min[1]=i;
	min_val = lyr[k].grid1d[i].cb();
      }
      if(i<lyr[k].numx()-1) i++;
      else {k++; i=0;}
    }
    if((real_min[0]==rmin[j][0] && real_min[1]==rmin[j][1]) ||(real_min[0]==lmin[j][0] && real_min[1]==lmin[j][1]) ) {
      for(l=j+1;l<n-1;l++) {
	min[l][0]=min[l+1][0];
	min[l][1]=min[l+1][1];
      }
      n--;
      j--;
    }

    j++;
  }

  /* Create dinamic arrays that will be returned */
  if(n>0) {
    left[0]  = new int [n];
    right[0] = new int [n];
    left[1]  = new int [n];
    right[1] = new int [n];
    for(j=0;j<n;j++) {
      left[0][j] =lmin[j][0];
      left[1][j] =lmin[j][1];
      right[0][j] =rmin[j][0];
      right[1][j] =rmin[j][1];
    }
  }
    
  return n;

}

		
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* FIND_HETERO_VB see FIND_HETERO_CB */

int Device::find_hetero_vb(int ** left, int ** right, FILE * fp) {

  const int MAX_STEPS = 10;
  const double band_step = 0.05;  /* eV/cm */

  int i,k,j,n,p,l,real_max[2],min[MAX_STEPS][2],lmin[MAX_STEPS][2],rmin[MAX_STEPS][2];
  double d,vb1,vb2,max_val;

  double max_length = sim->Max_quantum_width; /* cm */
  double secure_length = sim->Secure_length; /* cm */

  /* Find all the steps in the band structure */
  n = 0;
  vb1 = lyr[0].grid1d[0].chi()+lyr[0].grid1d[0].Eg();
  for(k=0;k<num_lyr;k++) {
    for(i=0;i<lyr[k].numx();i++) {
      vb2 = lyr[k].grid1d[i].chi()+lyr[k].grid1d[i].Eg();
      if( fabs(vb1-vb2)> band_step ) {
	min[n][0]=k;
	min[n][1]=i;
	n++;
      }
      vb1 = vb2;
    }
    if(n>=MAX_STEPS) {
      printf("FIND_HETERO_VB: Too many heterostructures\n");
      exit(1);
    }
  }

  if(n==0) {
    if(fp) printf("FIND_HETERO_VB: No heterostructures\n");
    return 0;
  }


  if(fp) printf("--------> %d steps\n",n);

  j = 0;
  /* LOOP ON STEPS */
  while(j<=n-1) {

    /* Set the left boundary position */
    if(min[j][1]==0) {
      lmin[j][0]= min[j][0]-1;
      lmin[j][1]= lyr[min[j][0]-1].numx()-1;
    }
    else {
      lmin[j][0]= min[j][0];
      lmin[j][1]= min[j][1]-1;
    }
    d = lyr[lmin[j][0]].grid1d[lmin[j][1]].xsize();
    while( d<secure_length && ( (j==0 && (lmin[j][0]>0 || (lmin[j][0]==0 && lmin[j][1]>0))) || (j>0 && ( lmin[j][0]>rmin[j-1][0] || ( lmin[j][0]==rmin[j-1][0] && lmin[j][1]>rmin[j-1][1] ))))) {
      if(lmin[j][1]==0) {
	lmin[j][0] = lmin[j][0]-1;
	lmin[j][1] = lyr[lmin[j][0]-1].numx()-1;
      }
      else lmin[j][1] = lmin[j][1]-1;
      d += lyr[lmin[j][0]].grid1d[lmin[j][1]].xsize();
    };
    if(lmin[j][0]==0 && lmin[j][1]==0) nrerror("FIND_HETERO_VB: The quantum structure has reached the left edge of the device. Decrease the secure_length variable or put an extra layer at the left of the device.");


    /* Check if within a max_length there is another step and change right boundary */
    rmin[j][0] = min[j][0];
    rmin[j][1] = min[j][1];
    k = min[j][0];
    i = min[j][1];
    d = 0.5*lyr[k].grid1d[i].xsize();
    if(i<lyr[k].numx()-1) i++;
    else {k++; i=0;}
    p = 1;
    while( d<max_length && k<num_lyr) {
      d += lyr[k].grid1d[i].xsize();
      if(i<lyr[k].numx()-1) i++;
      else {k++; i=0;}
      if(i==min[j+p][1] && k==min[j+p][0]) {
	rmin[j][0] = min[j+p][0];
	rmin[j][1] = min[j+p][1];
	d = 0;
	for(l=j+p;l<n-1;l++) {
	  min[l][0]=min[l+1][0];
	  min[l][1]=min[l+1][1];
	}
	n--;
      }
    }

    /* Set the right boundary position */
    d = 0;
    d += lyr[rmin[j][0]].grid1d[rmin[j][1]].xsize();
    while(d<secure_length && (rmin[j][0]<num_lyr-1 || ( rmin[j][0]==num_lyr-1 && rmin[j][1]<lyr[num_lyr-1].numx()-1))) {
      if(rmin[j][1]==lyr[rmin[j][0]].numx()-1) { 
	rmin[j][0] ++;
	rmin[j][1] = 0;
      }
      else rmin[j][1]++;
      d += lyr[rmin[j][0]].grid1d[rmin[j][1]].xsize();
    }
    if(rmin[j][0]==num_lyr-1 && rmin[j][1]==lyr[num_lyr-1].numx()-1)   nrerror("FIND_HETERO_VB: The quantum structure has reached the right edge of the device. Decrease the secure_length variable or put an extra layer at the right of the device.");


    /* Check whether there is a maximum */
    k = lmin[j][0];
    i = lmin[j][1];
    real_max[0]=k;
    real_max[1]=i;
    max_val = lyr[k].grid1d[i].vb();
    while(k<rmin[j][0] || (k==rmin[j][0] && i<=rmin[j][1])) {
      if(lyr[k].grid1d[i].vb()>max_val) {
	real_max[0]=k;
	real_max[1]=i;
	max_val = lyr[k].grid1d[i].vb();
      }
      if(i<lyr[k].numx()-1) i++;
      else {k++; i=0;}
    }
    if((real_max[0]==rmin[j][0] && real_max[1]==rmin[j][1]) ||(real_max[0]==lmin[j][0] && real_max[1]==lmin[j][1]) ) {
      for(l=j+1;l<n-1;l++) {
	min[l][0]=min[l+1][0];
	min[l][1]=min[l+1][1];
      }
      n--;
      j--;
    }

    j++;
  }


  /* Create dinamic arrays that will be returned */
  if(n>0) {
    left[0]  = new int [n];
    right[0] = new int [n];
    left[1]  = new int [n];
    right[1] = new int [n];
    for(j=0;j<n;j++) {
      left[0][j] =lmin[j][0];
      left[1][j] =lmin[j][1];
      right[0][j] =rmin[j][0];
      right[1][j] =rmin[j][1];
    }
  }
  return n;

}


/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* STEP_APPROXIMATION_CB creates an array of class "Step" called "cb_stp" and returns the number of 
   items insides this array. It uses left and right boundaries in the grid under the form (k=layer,i=position in the layer).
   The parameter "STEP_GAP" sets the minimal difference in energy between two steps.
   The parameter "MAX_STEP" is not necessary but it controls that the number of step is reasonably low.
   fp is the usual pointer to "quantum.log" file.
   fp1 points to a file in which the function dumps the step_approximation. If you don't desire to dump 
   any file, set the pointer to zero.  
class Step
   int           left[2];    position of the first(from left)mesh point in the step
   int           right[2];   position of the last(from left)mesh point in the step
   double        energy;     value in eV of the energy of the step.
   double        size;       lenght of the step in cm
   double        mass;       effective mass in the step
   Step *        ptr;        pointer to the next step 
*/ 
Step * Device::step_approximation_cb( int & num,int lf_lyr,int lf_pnt,int rh_lyr,int rh_pnt,FILE * fp,FILE * fp1 ) { 

  int i,k,l,n;
  double d;
  Step *top =0, *cr =0, * cb_stp, cr_tmp;

  const int MAX_STEP = 200;
  const double 	STEP_GAP = 1e-3;

  if(fp) fprintf(fp," STEP_APPROX.:CONDUCTION BAND  [%d,%d]<-->[%d,%d]\n",lf_lyr,lf_pnt,rh_lyr,rh_pnt);
  if(lf_lyr==rh_lyr && lf_pnt==rh_pnt) {
    printf("STEP_APPROX.:CONDUCTION BAND: No Steps\n");
    exit(1);
  }

  /* Create the chain of Step (top) */
  i   =	lf_pnt;
  k   =	lf_lyr;
  top = new Step;
  cr  =	top;
  cr->left[0] =	lf_lyr;
  cr->left[1] =	lf_pnt;
  cr->right[0] = lf_lyr;
  cr->right[1] = lf_pnt;
  cr->energy = lyr[k].grid1d[i].cb();
  cr->size = lyr[k].grid1d[i].xsize();
  cr->mass = lyr[k].grid1d[i].me();
  cr->ptr = 0;
  n = 1;

  while(k<rh_lyr || (i<rh_pnt && k==rh_lyr)) {
    if(i<lyr[k].numx()-1) i++;
    else { k++; i=0; }
    if(fabs(lyr[k].grid1d[i].cb()-(cr->energy))>STEP_GAP) {
      cr->ptr = new Step;
      if(cr->ptr == NULL) nrerror("Memory allocation failure in step_approximation_cb");
      cr = cr->ptr;
      cr->left[0] = k;
      cr->left[1] = i;
      cr->right[0] = k;
      cr->right[1] = i;
      cr->energy = lyr[k].grid1d[i].cb();
      cr->size = lyr[k].grid1d[i].xsize();
      cr->mass = lyr[k].grid1d[i].me();
      cr->ptr =	0;
      n++;
    }
    else {
      cr->right[0] = k;
      cr->right[1] = i;
      cr->size += lyr[k].grid1d[i].xsize();
    }
  }

  if(n>=MAX_STEP) {
    printf("STEP_APPROX.:CONDUCTION BAND:: Too many energy-steps in conduction band!\n");
    printf("1)Increase the value of MAX_STEP in main.h, or\n");
    printf("2)Increase the value of STEP_GAP in main.h!\n");
    exit(1);
  }

  /* Creates a dynamic vector(cb_stp) from the chain of Step(top) */
  cb_stp = new Step [n];
  if(cb_stp == NULL) nrerror("Memory allocation failure of cb_stp in step_approximation_cb");
  cr = top;
  for(i=0;i<n;i++) {
    cb_stp[i].left[0] = cr->left[0];
    cb_stp[i].left[1] = cr->left[1];
    cb_stp[i].right[0] = cr->right[0];
    cb_stp[i].right[1] = cr->right[1];
    cb_stp[i].energy = cr->energy;
    cb_stp[i].size = cr->size;
    cb_stp[i].mass = cr->mass;
    if(fp) fprintf(fp,"Step=%d, pos:(%d,%d)<--> (%d,%d) energy= %le, size= %le, mass= %lf\n"
	    ,i,cr->left[0],cr->left[1],cr->right[0],cr->right[1],cr->energy,cr->size,cr->mass);
    cr = cr->ptr;
  }


  /*Dump in the file the step approximation */
  if(fp1!=0) {
    d = 0;
    k=0;
    i=0;
    while(k!=lf_lyr || i!=lf_pnt) {
      d += lyr[k].grid1d[i].xsize();
      if(i<lyr[k].numx()-1) i++;
      else { k++; i=0; }
    }
    for(l=0;l<n;l++) {
      while(k<cb_stp[l].right[0] || ( k==cb_stp[l].right[0] && i<=cb_stp[l].right[1])) {
	d += lyr[k].grid1d[i].xsize();
	fprintf(fp1,"%le  %le\n",d,cb_stp[l].energy);
	if(i<lyr[k].numx()-1) i++;
	else { k++; i=0; }
      }
    }
  }

  /* Delete the chain of elements of type Step (top) */
  while(top!=NULL && top->ptr!=NULL) {
    cr = top;
    while(cr->ptr->ptr!=NULL) cr = cr->ptr;
    delete cr->ptr;
    cr->ptr = 0;
  }
  if(top!=NULL) delete top;

  num = n;
  return cb_stp;

}

/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */


Step * Device::step_approximation_vb( int & num,int lf_lyr,int lf_pnt,int rh_lyr,int rh_pnt,FILE * fp,FILE * fp1 ) { 

  int i,k,l,n;
  double d;
  Step * top, * cr, *vb_stp;

  const int MAX_STEP = 200;
  const double 	STEP_GAP = 1e-3;

  if(fp) fprintf( fp," STEP_APPROX.:VALENCE BAND  [%d,%d]<-->[%d,%d]\n",lf_lyr,lf_pnt,rh_lyr,rh_pnt);
  if(lf_lyr==rh_lyr && lf_pnt==rh_pnt) {
    printf("STEP_APPROX.:VALENCE BAND: No Steps\n");
    exit(1);
  }

  i = lf_pnt;
  k = lf_lyr;
  top = new Step;
  cr  =	top;
  cr->left[0] =	lf_lyr;
  cr->left[1] =	lf_pnt;
  cr->right[0] = lf_lyr;
  cr->right[1] = lf_pnt;
  cr->energy = -lyr[k].grid1d[i].vb();
  cr->size = lyr[k].grid1d[i].xsize();
  cr->mass = lyr[k].grid1d[i].mh();
  cr->ptr = 0;
  n = 1;

  while(k<rh_lyr || (i<rh_pnt && k==rh_lyr)) {
    if(i<lyr[k].numx()-1) i++;
    else { k++; i=0; }
    if(fabs(-lyr[k].grid1d[i].vb()-(cr->energy))>STEP_GAP) {
      cr->ptr =  new Step;
      if(cr->ptr == NULL) nrerror("Memory allocation failure in step_approximation_vb");
      cr = cr->ptr;
      cr->left[0] = k;
      cr->left[1] = i;
      cr->right[0] = k;
      cr->right[1] = i;
      cr->energy = -lyr[k].grid1d[i].vb();
      cr->size = lyr[k].grid1d[i].xsize();
      cr->mass = lyr[k].grid1d[i].mh();
      cr->ptr =	0;
      n++;
    }
    else {
      cr->right[0] = k;
      cr->right[1] = i;
      cr->size += lyr[k].grid1d[i].xsize();
    }
  }

  if(n>=MAX_STEP) {
    printf("STEP_APPROX.:VALENCE BAND:: Too many energy-steps in conduction band!\n");
    printf("1)Increase the value of MAX_STEP in main.h, or\n");
    printf("2)Increase the value of STEP_GAP in main.h!\n");
    exit(1);
  }

  vb_stp = new Step [n];
  if(vb_stp == NULL) nrerror("Memory allocation failure of vb_stp in step_approximation_vb");
  cr = top;
  for(i=0;i<n;i++) {
    vb_stp[i].left[0] 	= cr->left[0];
    vb_stp[i].left[1] 	= cr->left[1];
    vb_stp[i].right[0] 	= cr->right[0];
    vb_stp[i].right[1] 	= cr->right[1];
    vb_stp[i].energy 	= cr->energy;
    vb_stp[i].size 	= cr->size;
    vb_stp[i].mass 	= cr->mass;
    if(fp) fprintf(fp,"Step=%d, pos:(%d,%d)<--> (%d,%d) energy= %le, size= %le, mass= %lf\n"
	    ,i,cr->left[0],cr->left[1],cr->right[0],cr->right[1],cr->energy,cr->size,cr->mass);
    cr = cr->ptr;
  }

  /*Dump in the file the step approximation */
  if(fp1!=0) {
    d = 0;
    k=0;
    i=0;
    while(k!=lf_lyr || i!=lf_pnt) {
      d += lyr[k].grid1d[i].xsize();
      if(i<lyr[k].numx()-1) i++;
      else { k++; i=0; }
    }
    for(l=0;l<n;l++) {
      while(k<vb_stp[l].right[0] || ( k==vb_stp[l].right[0] && i<=vb_stp[l].right[1])) {
	d += lyr[k].grid1d[i].xsize();
	fprintf(fp1,"%le  %le\n",d,-vb_stp[l].energy);
	if(i<lyr[k].numx()-1) i++;
	else { k++; i=0; }
      }
    }
  }

  while(top!=NULL && top->ptr!=NULL) {
    cr = top;
    while(cr->ptr->ptr!=NULL) cr = cr->ptr;
    delete cr->ptr;
    cr->ptr = 0;
  }
  if(top!=NULL) delete top;
  num = n;
  return vb_stp;

}
//****************************************************************
//****************************************************************

/* Band is an array of num_step objects of type Step.
   E1 and E2 represents the initial and the final energy on which the transmission coefficient is calculated.
   gap is the step in energy from one to another evaluation of the transmission coefficient. 
   If gap<0 the transmission coefficient is calculated in MIN energy points from the band minimum energy untill the band maximum energy.
   If gap=0 the transmission coefficient is calculated in MIN energy points for OVER_TOP energy above the band maximum energy. 
   num_en contains the number of pairs (energy,transmission coefficient) returned by the function in a 2x1 array of doubles.
   fp is the usual pointer to the log file.
   left and right are the posistion of the initial step and final step of the array band on which the calculation is wanted to be done.
*/
double ** Device::transmission_coefficient(int num_step,Step * band,double E1,double E2,double gap,int & num_en,FILE * fp,int left,int right)	{

  int i,j,num;
  const int MIN = 500;
  const double OVER_TOP = 1.0; 
  double d,E_max,E_min,E,E_fin,E_in,E_gap;

  double ** tr = 0;
  Complex * k_vec = 0,prod;
  Tmatrix * Tmat = 0;


  if(right>left) num = right - left +1;
  else num = num_step;

  if(fp) fprintf(fp,"Calculation of the transmission coefficient ...\n");

  E_gap = gap;
  E_in 	= E1;
  E_fin = E2;

  k_vec	= new Complex [num];
  if(k_vec==NULL) nrerror("Memory allocation failure for k_vec in transmission_coefficient");
  Tmat = new Tmatrix [num-1];
  if(Tmat==NULL) nrerror("Memory allocation failure for Tmat in transmission_coefficient");

  /* SET THE TYPE OF CALCULATION */
  if(E_gap<0) {
    E_max = -50.0;
    E_min = 50.0;
    for(j=0;j<num;j++) {
      if(E_max<band[j+left].energy)	E_max=band[j+left].energy;
      if(E_min>band[j+left].energy)	E_min=band[j+left].energy;
    }
    E_in  = E_min;
    E_fin = E_max;
    num_en= MIN;
    E_gap = (E_fin - E_in )/(double)MIN;
  }
  else if(E_gap==0) {
    E_max = -50.0;
    for(j=0;j<num;j++) if(E_max<band[j+left].energy) E_max=band[j+left].energy;
    E_in  = E_max;
    E_fin = E_max+OVER_TOP;
    num_en= MIN;
    E_gap = OVER_TOP/(double)MIN;
  }
  else if(E_fin==E_in) num_en = 1;
  else num_en = int((double)(E_fin-E_in )/E_gap);

  /* Calculation at least over MIN points */
  if(num_en<MIN) {
    num_en = MIN;
    E_gap  = (E_fin - E_in )/(double)MIN;
  }

  /* ALLOCATION ARRAYS */
  tr = new double* [2];
  if(tr==NULL) nrerror("Memory allocation failure for tr in transmission_coefficient");
  tr[0] = new double [num_en];
  if(tr[0]==NULL) nrerror("Memory allocation failure for tr[0] in transmission_coefficient");
  tr[1] = new double [num_en];
  if(tr[1]==NULL) nrerror("Memory allocation failure for tr[1] in transmission_coefficient");

  E = E_in + 0.1*E_gap;
  for(i=0;i<num_en;i++)	{ // Iteration on Energy
    for(j=0;j<num;j++) { // Iteration on step.
      if(fabs(E-band[j+left].energy)<1.0e-20) {
	printf("Transmission_coefficient: Energy equals to step: correction made!\n");
	E = E+0.01*E_gap;
      }
      k_vec[j]=cplx_sqrt(band[j+left].mass*KKK*(E-band[j+left].energy));
    }
    d =	0.0;
    for(j=0;j<num-1;j++)	{
      Tmat[j].is(k_vec[j],k_vec[j+1],band[j+left].mass,band[j+1+left].mass,d);
      d	+= band[j+1+left].size*LNTHCONV;
    }
			
    prod = product(num-1,Tmat);

    if(prod.rl==0 && prod.im==0) {
      /* printf("CALC_TRANSMISSION: Zero transmission coefficient at %d  E=%le in [%le,%le]\n",i,E, E_in,E_fin);
	 delete [] tr[0];
	 delete [] tr[1];
	 delete [] tr;
	 delete [] k_vec;
	 delete [] Tmat;
	 exit(1); 
      */
      if(i>0)  tr[1][i]=tr[1][i-1];
      else tr[1][i]=1e-20;
    }
    else tr[1][i]=cplx_mod_sqr(prod);

    tr[0][i]=E;
    E += E_gap;

  }

  delete [] k_vec;
  delete [] Tmat;

  if(fp) fprintf(fp,"Transmission coefficient done!\n");

  return tr;

}



/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */


void Device::wavefunction_n( Eigenvalue ** eigen_store,int step_num,Step * band,FILE * fp, char str1[30]) {

  short switch1 = 1; // Check the correct curvature of the wave function.
  short switch2 = 1; // Check wave functions with similar eigenvalue.
  short switch3 = 1; // Chech for non-negligble tails of wavefunctions
  short switch4 = 1; // Check for the lowest energy value.

  int i,k,j,j_pos,t,ll,k_old,i_old,count,ndone,ckr,Ebuffer[10];
  int nsim,min_prob,bug_up,peak_count,found,n_max;
  double d,prob1,E,bb,sum,d0,done[50],right_prob[10],prob2,prob3,h,E_min,E_max,old_sum;
  const double  ENERGY_TOLLERANCE = 0.0005; // Max distance of similar eigenvalues (eV)
  const double  PROB_TAIL = 1e-3; // Max probability accepted in the tail of the wavefunction
  Complex	Ac,Bc;
  Qntm_var **    ptr_qntm;
  Eigenvalue *  tmp;
  FILE *        fp1;

  double secure_length = sim->Secure_length; /* cm */

  if(fp) printf("WAVEFUNCTION_N STARTING ...\n");
  if(step_num == 0) nrerror("Wavefunction_n: num of steps = 0");

  if(fp) printf("It calculates the distance to which we have the heterostructure\n");
  k = 0;
  i = 0;
  d = 0;
  while(test_upper_limit(k,i,band[0].left[0],band[0].left[1])) {
    d += lyr[k].grid1d[i].xsize();
    lyr[k].grid1d[i].qptr_e = 0;
    incr(k,i);
  }
  d0 = d;
  decr(k,i);
  d -= lyr[k].grid1d[i].xsize();

  if(fp) printf("It creates and initilizes every Qntm_var class for each mesh point of the heterostructure \n"); 
  j = 0;
  while(j<step_num) {
    if(k!=band[j].left[0] || i != band[j].left[1]) nrerror("Wavefunction_n:: error in mesh scanning");
    while(test_upper_limit(k,i,band[j].right[0],band[j].right[1]) ) {
      if(lyr[k].grid1d[i].qptr_e!=NULL)  nrerror("Wavefunction_n: writing on a previous quantum structure");
      lyr[k].grid1d[i].qptr_e = new Qntm_var;
      lyr[k].grid1d[i].qptr_e->set_eigen(eigen_store[j],band[j].energy,d,d0,lyr[k].grid1d[i].me());
      d += lyr[k].grid1d[i].xsize();
      incr(k,i);
    }
    j++;
  }
  while(test_upper_limit(k,i)) {
    lyr[k].grid1d[i].qptr_e = 0;
    incr(k,i);
  }


  if(fp) printf("Normalization ... \n");
  ll = 0;
  ndone = 0; /* Number of wavefunction already calculate of energy given by done[] */
  while(ll<step_num) { /* loop on steps */
    tmp	= eigen_store[ll];
    while(tmp!=0) {
      E	= tmp->energy;

      /* check if that eigenvalue energy has been calculated before */
      ckr = ndone-1;
      while(ckr>=0 && done[ckr]!=E) ckr--;
      if(ckr<0) {
	done[ndone]=E;
	ndone++;
	if(fp) printf("Normalization of the wave function\n");
	if(!electron_wavefunction_normalization(E,fp)) nrerror("Wavefunction_n:: Not normalized!");
	
	//-------------------------------------------------------
	if(switch1) { 
	  if(fp) printf("Check the correct curvature\n");
	  do{
	    bug_up = 0;
	    k=0;
	    i=0;
	    while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e==0 ) incr(k,i);
	    while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->check_norm(E)!=-1 ) incr(k,i);
	    decr(k,i);
	    k_old = k;
	    i_old = i;
	    decr(k,i);
	    while(lyr[k].grid1d[i].qptr_e && lyr[k_old].grid1d[i_old].qptr_e &&
		  (lyr[k].grid1d[i].qptr_e->calc_prob(E,0)<lyr[k_old].grid1d[i_old].qptr_e->calc_prob(E,0))) {
	      if(fp) printf("WAVE_FUNCTION_N: Tail going up (%d,%d)!\n",k,i);
	      bug_up = 1;
	      if(!lyr[k_old].grid1d[i_old].qptr_e->drop(E)) nrerror("Wavefunction_n: cannot drop eigenvalue");
	      k_old = k;
	      i_old = i;
	      decr(k,i);
	    }
	    if(fp) printf("Normalization of the wave function\n");
	    if(!electron_wavefunction_normalization(E,fp)) nrerror("Wavefunction_n:: Not normalized!");
	  } while(bug_up);
	}
	// -------------------------------------------------------

      }
      tmp = tmp->ptr;
    }
    ll++;
  }


  if(fp) printf("Ordering in crescent values of energies\n");
  ckr = 1;
  while(ckr) {
    ckr = 0;
    for(ll=1;ll<ndone;ll++) {
      if(done[ll-1]>done[ll]) {
	E = done[ll];
	done[ll] = done[ll-1];
	done[ll-1] = E;
	ckr = 1;
      }
    }
  }

  // ---------------------------------------------------------------
  if(switch2) {
    if(fp) printf("Check whether there are similar wavefunctions\n");
    nsim = 0;
    for(j=1;j<ndone;j++) {
      if(done[j]-done[j-1]<ENERGY_TOLLERANCE) Ebuffer[nsim++]=j-1;
      else {
	if(nsim>0) {
	  if(fp) printf("Similar wavefunctions found!\n");
	  Ebuffer[nsim++]=j-1;
	  min_prob = 0;
	  for(ll=0;ll<nsim;ll++) {
	    E = done[Ebuffer[ll]];
	    k = band[0].left[0];
	    i = band[0].left[1];
	    while(lyr[k].grid1d[i].qptr_e->calc_prob(E,0)==0) {
	      if(i<lyr[k].numx()-1) i++;
	      else { k++; i = 0; }
	    }
	    while( (k<band[step_num-1].right[0] || (k==band[step_num-1].right[0] && i<band[step_num-1].right[1])) 
		   && lyr[k].grid1d[i].qptr_e->calc_prob(E,0)>0 ) {
	      k_old = k;
	      i_old = i;
	      if(i<lyr[k].numx()-1) i++;
	      else { k++; i = 0; }
	    }
	    right_prob[ll] = lyr[k_old].grid1d[i_old].qptr_e->calc_prob(E,0);
	    if(right_prob[ll]<right_prob[min_prob]) min_prob = ll; 
	  }
	  for(ll=0;ll<nsim;ll++) 
	    if(ll!=min_prob) {
	      k = band[0].left[0];
	      i = band[0].left[1];
	      while( k<band[step_num-1].right[0] || (k==band[step_num-1].right[0] && i<=band[step_num-1].right[1]) ) {
		if(lyr[k].grid1d[i].qptr_e->drop(done[Ebuffer[ll]])) {
		  delete lyr[k].grid1d[i].qptr_e;
		  lyr[k].grid1d[i].qptr_e = 0;
		}
		if(i<lyr[k].numx()-1) i++;
		else { k++; i = 0; }
	      }
	    }
	  count = 0;
	  for(ll=0;ll<nsim;ll++) {
	    if(ll!=min_prob) {
	      for(t=Ebuffer[ll]-count;t<ndone;t++) done[t]=done[t+1];
	      ndone--;
	      count++;
	    }
	  }
	  nsim = 0;
	}
      }
    }
  }
  // --------------------------------------------------------------


  if(switch3) {
    if(fp) printf("Check non-negligible wavefunction tails\n");
    for(j=0;j<ndone;j++) {
      
      E = done[j];
      if(fp) printf("Left tail at E=%lf...\n",E);
      k = 0;
      i = 0;
      while(test_upper_limit(k,i) && (!(lyr[k].grid1d[i].qptr_e) || lyr[k].grid1d[i].qptr_e->calc_prob(E,0)==0)) incr(k,i);
      if(test_upper_limit(k,i) && (lyr[k].grid1d[i].qptr_e)!=0) prob1 = lyr[k].grid1d[i].qptr_e->calc_prob(E,0);
      else prob1 = 0;
      if(prob1>PROB_TAIL) {
	d = 0;
	while(d<secure_length && test_lower_limit(k,i)) {
	  k_old = k;
	  i_old = i;
	  decr(k,i);
	  d += lyr[k].grid1d[i].xsize();
	  if(lyr[k].grid1d[i].qptr_e==0) {
	    if(fp) printf("Adding new Qntm_var structure in (%d,%d)\n",k,i);
	    lyr[k].grid1d[i].qptr_e = new Qntm_var;
	    (* lyr[k].grid1d[i].qptr_e) = (* lyr[k_old].grid1d[i_old].qptr_e);
	    lyr[k].grid1d[i].qptr_e ->set_offset(lyr[k_old].grid1d[i_old].qptr_e->ret_offset()-lyr[k].grid1d[i].xsize()*LNTHCONV);
	  }
	  else {
	    if(fp) printf("Adding an eigenvalue to the pre-existed structure\n");
	    if(!lyr[k].grid1d[i].qptr_e->add_eigenvalue( lyr[k_old].grid1d[i_old].qptr_e->ret_eval(E),
							 lyr[k_old].grid1d[i_old].qptr_e->ret_wave(E) )) {
	      lyr[k].grid1d[i].qptr_e->modify_wave( lyr[k_old].grid1d[i_old].qptr_e->ret_eval(E),
						    lyr[k_old].grid1d[i_old].qptr_e->ret_wave(E) );
	      printf("WAVEFUNCTION_N:1 Eigenvalue already present\n");
	    }
	  }
	  if(lyr[k].grid1d[i].qptr_e->calc_prob(E,0)>lyr[k_old].grid1d[i_old].qptr_e->calc_prob(E,0)) d = secure_length ;
	}
      }
      

      if(fp) printf("Right tail at E=%lf...\n",E);
      k = 0;
      i = 0;
      while(k<=band[step_num-1].right[0] && (!(lyr[k].grid1d[i].qptr_e) || lyr[k].grid1d[i].qptr_e->calc_prob(E,0)==0)) incr(k,i);
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e!=0  && lyr[k].grid1d[i].qptr_e->calc_prob(E,0)>0) incr(k,i);
      if(test_upper_limit(k,i)) {
	k_old = k;
	i_old = i;
	decr(k_old,i_old);
	prob1 = lyr[k_old].grid1d[i_old].qptr_e->calc_prob(E,0);
	if(prob1>PROB_TAIL) {
	  d = lyr[k].grid1d[i].xsize();
	  while(d<secure_length && test_upper_limit(k,i) && prob1>PROB_TAIL) {
	    if(lyr[k].grid1d[i].qptr_e==0) {
	      if(fp) printf("Adding new Qntm_var structure in (%d,%d)\n",k,i);
	      lyr[k].grid1d[i].qptr_e = new Qntm_var;
	      if(!lyr[k].grid1d[i].qptr_e) nrerror("WAVEFUNCTION_N: Memory allocation failure!");
	      ( * lyr[k].grid1d[i].qptr_e) =  (* lyr[k_old].grid1d[i_old].qptr_e);
	      lyr[k].grid1d[i].qptr_e ->set_offset(lyr[k_old].grid1d[i_old].qptr_e->ret_offset()+lyr[k].grid1d[i].xsize()*LNTHCONV);
	    }
	    else {
	      if(fp) printf("Adding an eigenvalue to the pre-existed structure\n");
	      if(!lyr[k].grid1d[i].qptr_e->add_eigenvalue( lyr[k_old].grid1d[i_old].qptr_e->ret_eval(E),
							   lyr[k_old].grid1d[i_old].qptr_e->ret_wave(E) )) {
		lyr[k].grid1d[i].qptr_e->modify_wave( lyr[k_old].grid1d[i_old].qptr_e->ret_eval(E),
						      lyr[k_old].grid1d[i_old].qptr_e->ret_wave(E) );
		if(fp) printf("WAVEFUNCTION_N:2 Eigenvalue already present\n");
	      }
	    }
	    if(lyr[k].grid1d[i].qptr_e->calc_prob(E,0)>lyr[k_old].grid1d[i_old].qptr_e->calc_prob(E,0)) d = secure_length ;
	    k_old = k;
	    i_old = i;
	    incr(k,i);
	    d += lyr[k].grid1d[i].xsize();
	  }
	}
      }
      
      if(!electron_wavefunction_normalization(E,fp)) nrerror("Wavefunction_n:: Not normalized!");
    }
  }



  // ----------------------------------------------------------------------
  if(switch4) {
    if(fp) printf("Find the lowest energy wavefunction\n");
    k = 0;
    i = 0;
    E = 1e10;
    while(test_upper_limit(k,i)) {
      if(lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->ret_num()>0 && lyr[k].grid1d[i].qptr_e->ret_energy(0)<E )
	E = lyr[k].grid1d[i].qptr_e->ret_energy(0);
      incr(k,i);
    }

    if(E<1e10) {
      if(fp) printf("Spot the number of peaks of the wave function at energy %lf\n",E);
      peak_count = 0;
      j = 0;
      k = 0;
      i = 0;
      prob1 = -1.0;
      prob2 = -1.0;
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e==0 ) incr(k,i);
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->check_norm(E)==-1 ) incr(k,i);
      if(lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->check_norm(E)!=-1) prob1=lyr[k].grid1d[i].qptr_e->calc_prob(E,0);
      incr(k,i);
      if(lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->check_norm(E)!=-1) prob2=lyr[k].grid1d[i].qptr_e->calc_prob(E,0);
      incr(k,i);
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->check_norm(E)!=-1 && prob1!=-1.0 && prob2!=-1.0) {
	prob3=lyr[k].grid1d[i].qptr_e->calc_prob(E,0);
	if(prob1<prob2 && prob2>prob3) peak_count++;
	prob1=prob2;
	prob2=prob3;
	incr(k,i);
      }
      if(fp) printf(" .... %d peaks found\n",peak_count);
    }
    else {
      if(fp) printf("It seems there is not a single wavefunction.\n");
      peak_count = 2;      
    }
    
    if(fp) printf("Find the bottom of conduction band\n");
    j = 1;
    j_pos = 0;
    E_min = 999.0;
    while(j<step_num-1) {
      if(band[j].energy<band[j-1].energy && band[j].energy<band[j+1].energy && band[j].energy<E_min) {
	E_min = band[j].energy;
	j_pos = j;
      }
      j++;
    }

    if(peak_count>1 && j_pos>0) {
      
      if(E==1e10) {
	j = 0;
	while(band[j_pos-j].energy<band[j_pos-j-1].energy && band[j_pos+j].energy<band[j_pos+j+1].energy) j++;
	E_max = (band[j_pos-j].energy>band[j_pos+j+1].energy)?band[j_pos+j].energy:band[j_pos-j].energy;
      }
      else E_max = E;
      if(E_min>=E_max) nrerror("WAVEFUNCTION_N1: first wavefunction not found!");
      
      h = (E_max - E_min)*1e-4;
      E = E_min+h;
      found = 0;
      old_sum = 1e250;

      while(E<E_max && found !=1 ) {

	if(fp) printf("It creates and initilizes every Qntm_var class for each mesh point of the heterostructure\n"); 
	ptr_qntm = new Qntm_var * [xnum_mesh_points];
	tmp = new Eigenvalue;
	tmp->energy = E;
	tmp->gamma = 0;
	tmp->sigma = 0;
	tmp->ptr = 0;
	t = 0;
	k = 0;
	i = 0;
	while(test_upper_limit(k,i)) {
	  if(lyr[k].grid1d[i].qptr_e!=0 ) { // && lyr[k].grid1d[i].qptr_e->check_norm(E)!=-1
	    ptr_qntm[t] = new Qntm_var;
	    ptr_qntm[t]->set_eigen(tmp, lyr[k].grid1d[i].qptr_e->ret_Eo(),
				   (lyr[k].grid1d[i].qptr_e->ret_offset()+lyr[k].grid1d[i].qptr_e->ret_leftedge())/LNTHCONV,
				   lyr[k].grid1d[i].qptr_e->ret_leftedge()/LNTHCONV,lyr[k].grid1d[i].me());
	  }
	  else ptr_qntm[t] = 0;
	  incr(k,i);
	  t++;
	}

	if(fp) printf("Set the complex coefficients of the wave function\n");
	bb = 1.0;
	sum = 0.0;
	k = 0;
	i = 0;
	t = 0;
	while(test_upper_limit(k,i) && ptr_qntm[t]==0 ) {
	  incr(k,i);
	  t++;
	}
	if(test_upper_limit(k,i)) {
	  Ac.rl = 0.0;
	  Ac.im = 0.0;
	  Bc.rl = bb;
	  Bc.im = 0.0;
	  if(!ptr_qntm[t]->set_wave(E,Ac,Bc)) nrerror("N_NORMALIZTION out of sync");
	  sum += ptr_qntm[t]->calc_prob(E,lyr[k].grid1d[i].xsize());
	  incr(k,i);
	  t++;
	  while( ptr_qntm[t]!=0 && ptr_qntm[t]->check_norm(E)!=-1) {
	    if(fabs(ptr_qntm[t]->ret_Eo()-ptr_qntm[t-1]->ret_Eo())>1e-6)
	      ptr_qntm[t]->propagation(E,Ac,Bc,ptr_qntm[t-1]);
	    ptr_qntm[t]->set_wave(E,Ac,Bc);
	    sum += ptr_qntm[t]->calc_prob(E,lyr[k].grid1d[i].xsize());
	    incr(k,i);
	    t++;
	  }
	}
	
	if(fp) printf("Spotting the number of peaks\n");
	n_max = 0;
	t = 1;
	while(t<xnum_mesh_points-1) {
	  ckr = 0;
	  if(ptr_qntm[t-1] && ptr_qntm[t] && ptr_qntm[t+1]) {
	    if(ptr_qntm[t]->calc_prob(E,0)>ptr_qntm[t-1]->calc_prob(E,0)
	       && ptr_qntm[t]->calc_prob(E,0)>ptr_qntm[t+1]->calc_prob(E,0)) ckr = 1;
	  }
	  if(ckr) n_max++;
	  t++;
	}
	if(fp) printf("In E=%lf I found %d peaks\n",E,n_max);

	if(n_max==1 && sum>old_sum) {
	  fopen_s(&fp1,"OUTPUT_FILES/tempo_n.dat","w");
	  t = 0;
	  while(t<xnum_mesh_points) {
	    if(ptr_qntm[t]) fprintf(fp1,"%d  %le\n",t,ptr_qntm[t]->calc_prob(E,0));
	    t++;
	  }
	  fclose(fp1);
	  found = 1;
	  E -= h;
	}
	else {
	  old_sum = sum;
	  E += h;
	}

	if(fp) printf("Removing temporary arrays\n"); 
	delete tmp;
	for(i=0;i<xnum_mesh_points;i++) {
	  if(ptr_qntm[i]) delete ptr_qntm[i];
	  ptr_qntm[i] = 0;
	}
	delete [] ptr_qntm;
	ptr_qntm = 0;
	
      }
      
      if(found) {

	if(fp) printf("Adding new eigenvalue\n"); 
	j = 0;
	k = 0;
	i = 0;
	t = 0;
	while(test_upper_limit(k,i)) {
	  if(lyr[k].grid1d[i].qptr_e) lyr[k].grid1d[i].qptr_e->add_eigenvalue(E,lyr[k].grid1d[i].me());
	  incr(k,i);
	}
	if(fp) printf("Normalization of the wave function\n");
	if(!electron_wavefunction_normalization(E,fp)) nrerror("Wavefunction_n:: Not normalized!");
	
	for(ll=ndone-1;ll>=0;ll--) done[ll+1]=done[ll];
	done[0]=E;
	ndone++;
      }
      else printf("WAVEFUNCTION_N: Lowest value not found!\n");
      
    }
  }



  if(fp) printf("Set the probability vector and dump it in a file\n");
  for(ll=0;ll<ndone;ll++) {    
    E = done[ll];
    fopen_s(&fp1,str1,"w");
    j = 0;
    k = band[0].left[0];
    i = band[0].left[1];
    while(j<step_num) {
      while(test_upper_limit(k,i,band[j].right[0],band[j].right[1]) ) {
	lyr[k].grid1d[i].qptr_e->set_prob_e(E,lyr[k].grid1d[i].xsize(),fp1);
	incr(k,i);
      }
      j++;
    }
    fclose(fp1);
    str1[19]++;
  }


}





/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */


void Device::wavefunction_p( Eigenvalue ** eigen_store,int step_num,Step * band,FILE * fp, char str1[30]) {

  short switch1 = 1; // Check the correct curvature of the wave function.
  short switch2 = 1; // Check wave functions with similar eigenvalue.
  short switch3 = 1; // Chech for non-negligble tails of wavefunctions
  short switch4 = 1; // Check for the lowest energy value.

  int 		i,k,j,j_pos,t,ll,k_old,i_old,count,ndone,ckr,Ebuffer[10];
  int           nsim,min_prob,bug_up,peak_count,found,n_max;
  double 	d,prob1,E,bb,sum,d0,done[50],right_prob[10],prob2,prob3,h,E_min,E_max,old_sum;
  const double  ENERGY_TOLLERANCE = 0.0005; // Max distance of similar eigenvalues (eV)
  const double  PROB_TAIL = 1e-3; // Max probability accepted in the tail of the wavefunction
  Complex	Ac,Bc;
  Qntm_var **    ptr_qntm;
  Eigenvalue *  tmp;
  FILE *        fp1;

  double secure_length = sim->Secure_length; /* cm */

  if(fp) printf("WAVEFUNCTION_P STARTING ...\n");
  if(step_num == 0) nrerror("Wavefunction_p: num of steps = 0");

  if(fp) printf("It calculates the distance to which we have the heterostructure\n");
  k = 0;
  i = 0;
  d = 0;
  while(test_upper_limit(k,i,band[0].left[0],band[0].left[1])) {
    d += lyr[k].grid1d[i].xsize();
    lyr[k].grid1d[i].qptr_h = 0;
    incr(k,i);
  }
  d0 = d;
  decr(k,i);
  d -= lyr[k].grid1d[i].xsize();

  if(fp) printf("It creates and initilizes every Qntm_var class for each mesh point of the heterostructure \n"); 
  j = 0;
  while(j<step_num) {
    if(k!=band[j].left[0] || i != band[j].left[1]) nrerror("Wavefunction_p:: error in mesh scanning");
    while(test_upper_limit(k,i,band[j].right[0],band[j].right[1]) ) {
      if(lyr[k].grid1d[i].qptr_h!=NULL)  nrerror("Wavefunction_p: writing on a previous quantum structure");
      lyr[k].grid1d[i].qptr_h = new Qntm_var;
      lyr[k].grid1d[i].qptr_h->set_eigen(eigen_store[j],band[j].energy,d,d0,lyr[k].grid1d[i].mh());
      d += lyr[k].grid1d[i].xsize();
      incr(k,i);
    }
    j++;
  }
  while(test_upper_limit(k,i)) {
    lyr[k].grid1d[i].qptr_h = 0;
    incr(k,i);
  }


  if(fp) printf("Normalization ... \n");
  ll = 0;
  ndone = 0; /* Number of wavefunction already calculate of energy given by done[] */
  while(ll<step_num) { /* loop on steps */
    tmp	= eigen_store[ll];
    while(tmp!=0) {
      E	= tmp->energy;

      /* check if that eigenvalue energy has been calculated before */
      ckr = ndone-1;
      while(ckr>=0 && done[ckr]!=E) ckr--;
      if(ckr<0) {
	done[ndone]=E;
	ndone++;
	if(fp) printf("Normalization of the wave function\n");
	if(!hole_wavefunction_normalization(E,fp)) nrerror("Wavefunction_p:: Not normalized!");
	
	//-------------------------------------------------------
	if(switch1) { 
	  if(fp) printf("Check the correct curvature\n");
	  do{
	    bug_up = 0;
	    k = 0;
	    i = 0;
	    while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e==0 ) incr(k,i);
	    while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->check_norm(E)!=-1 ) incr(k,i);
	    decr(k,i);
	    k_old = k;
	    i_old = i;
	    decr(k,i);
	    while(lyr[k].grid1d[i].qptr_h && lyr[k_old].grid1d[i_old].qptr_h &&
		  (lyr[k].grid1d[i].qptr_h->calc_prob(E,0)<lyr[k_old].grid1d[i_old].qptr_h->calc_prob(E,0))) {
	      if(fp) printf("WAVE_FUNCTION_P: Tail going up (%d,%d)!\n",k,i);
	      bug_up = 1;
	      if(!lyr[k_old].grid1d[i_old].qptr_h->drop(E)) nrerror("Wavefunction_p: cannot drop eigenvalue");
	      k_old = k;
	      i_old = i;
	      decr(k,i);
	    }
	    if(fp) printf("Normalization of the wave function\n");
	    if(!hole_wavefunction_normalization(E,fp)) nrerror("Wavefunction_p:: Not normalized!");
	  } while(bug_up);
	}
	// -------------------------------------------------------

      }
      tmp = tmp->ptr;
    }
    ll++;
  }


  if(fp) printf("Ordering in crescent values of energies\n");
  ckr = 1;
  while(ckr) {
    ckr = 0;
    for(ll=1;ll<ndone;ll++) {
      if(done[ll-1]>done[ll]) {
	E = done[ll];
	done[ll] = done[ll-1];
	done[ll-1] = E;
	ckr = 1;
      }
    }
  }

  // ---------------------------------------------------------------
  if(switch2) {
    if(fp) printf("Check whether there are similar wavefunctions\n");
    nsim = 0;
    for(j=1;j<ndone;j++) {
      if(done[j]-done[j-1]<ENERGY_TOLLERANCE) Ebuffer[nsim++]=j-1;
      else {
	if(nsim>0) {
	  if(fp) printf("Similar wavefunctions found!\n");
	  Ebuffer[nsim++]=j-1;
	  min_prob = 0;
	  for(ll=0;ll<nsim;ll++) {
	    E = done[Ebuffer[ll]];
	    k = band[0].left[0];
	    i = band[0].left[1];
	    while(lyr[k].grid1d[i].qptr_h->calc_prob(E,0)==0) {
	      if(i<lyr[k].numx()-1) i++;
	      else { k++; i = 0; }
	    }
	    while( (k<band[step_num-1].right[0] || (k==band[step_num-1].right[0] && i<band[step_num-1].right[1])) 
		   && lyr[k].grid1d[i].qptr_h->calc_prob(E,0)>0 ) {
	      k_old = k;
	      i_old = i;
	      if(i<lyr[k].numx()-1) i++;
	      else { k++; i = 0; }
	    }
	    right_prob[ll] = lyr[k_old].grid1d[i_old].qptr_h->calc_prob(E,0);
	    if(right_prob[ll]<right_prob[min_prob]) min_prob = ll; 
	  }
	  for(ll=0;ll<nsim;ll++) 
	    if(ll!=min_prob) {
	      k = band[0].left[0];
	      i = band[0].left[1];
	      while( k<band[step_num-1].right[0] || (k==band[step_num-1].right[0] && i<=band[step_num-1].right[1]) ) {
		if(lyr[k].grid1d[i].qptr_h->drop(done[Ebuffer[ll]])) {
		  delete lyr[k].grid1d[i].qptr_h;
		  lyr[k].grid1d[i].qptr_h = 0;
		}
		if(i<lyr[k].numx()-1) i++;
		else { k++; i = 0; }
	      }
	    }
	  count = 0;
	  for(ll=0;ll<nsim;ll++) {
	    if(ll!=min_prob) {
	      for(t=Ebuffer[ll]-count;t<ndone;t++) done[t]=done[t+1];
	      ndone--;
	      count++;
	    }
	  }
	  nsim = 0;
	}
      }
    }
  }
  // --------------------------------------------------------------


  if(switch3) {
    if(fp) printf("Check non-negligible wavefunction tails\n");
    for(j=0;j<ndone;j++) {
      
      E = done[j];
      if(fp) printf("Left tail at E=%lf...\n",E);
      k = 0;
      i = 0;
      while(test_upper_limit(k,i) && (!(lyr[k].grid1d[i].qptr_h) || lyr[k].grid1d[i].qptr_h->calc_prob(E,0)==0)) incr(k,i);
      if(test_upper_limit(k,i) && (lyr[k].grid1d[i].qptr_h)!=0) prob1 = lyr[k].grid1d[i].qptr_h->calc_prob(E,0);
      else prob1 = 0;
      if(prob1>PROB_TAIL) {
	d = 0;
	while(d<secure_length && test_lower_limit(k,i)) {
	  k_old = k;
	  i_old = i;
	  decr(k,i);
	  d += lyr[k].grid1d[i].xsize();
	  if(lyr[k].grid1d[i].qptr_h==0) {
	    if(fp) printf("Adding new Qntm_var structure in (%d,%d)\n",k,i);
	    lyr[k].grid1d[i].qptr_h = new Qntm_var;
	    (* lyr[k].grid1d[i].qptr_h) = (* lyr[k_old].grid1d[i_old].qptr_h);
	    lyr[k].grid1d[i].qptr_h ->set_offset(lyr[k_old].grid1d[i_old].qptr_h->ret_offset()-lyr[k].grid1d[i].xsize()*LNTHCONV);
	  }
	  else {
	    if(fp) printf("Adding an eigenvalue to the pre-existed structure\n");
	    if(!lyr[k].grid1d[i].qptr_h->add_eigenvalue( lyr[k_old].grid1d[i_old].qptr_h->ret_eval(E),
							 lyr[k_old].grid1d[i_old].qptr_h->ret_wave(E) )) {
	      lyr[k].grid1d[i].qptr_h->modify_wave( lyr[k_old].grid1d[i_old].qptr_h->ret_eval(E),
						    lyr[k_old].grid1d[i_old].qptr_h->ret_wave(E) );
	      printf("WAVEFUNCTION_P:1 Eigenvalue already present\n");
	    }
	  }
	  if(lyr[k].grid1d[i].qptr_h->calc_prob(E,0)>lyr[k_old].grid1d[i_old].qptr_h->calc_prob(E,0)) d = secure_length ;
	}
      }
      

      if(fp) printf("Right tail at E=%lf...\n",E);
      k = 0;
      i = 0;
      while(k<=band[step_num-1].right[0] && (!(lyr[k].grid1d[i].qptr_h) || lyr[k].grid1d[i].qptr_h->calc_prob(E,0)==0)) incr(k,i);
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_h!=0  && lyr[k].grid1d[i].qptr_h->calc_prob(E,0)>0) incr(k,i);
      if(test_upper_limit(k,i)) {
	k_old = k;
	i_old = i;
	decr(k_old,i_old);
	prob1 = lyr[k_old].grid1d[i_old].qptr_h->calc_prob(E,0);
	if(prob1>PROB_TAIL) {
	  d = lyr[k].grid1d[i].xsize();
	  while(d<secure_length && test_upper_limit(k,i) && prob1>PROB_TAIL) {
	    if(lyr[k].grid1d[i].qptr_h==0) {
	      if(fp) printf("Adding new Qntm_var structure in (%d,%d)\n",k,i);
	      lyr[k].grid1d[i].qptr_h = new Qntm_var;
	      if(!lyr[k].grid1d[i].qptr_h) nrerror("WAVEFUNCTION_P: Memory allocation failure!");
	      ( * lyr[k].grid1d[i].qptr_h) =  (* lyr[k_old].grid1d[i_old].qptr_h);
	      lyr[k].grid1d[i].qptr_h ->set_offset(lyr[k_old].grid1d[i_old].qptr_h->ret_offset()+lyr[k].grid1d[i].xsize()*LNTHCONV);
	    }
	    else {
	      if(fp) printf("Adding an eigenvalue to the pre-existed structure\n");
	      if(!lyr[k].grid1d[i].qptr_h->add_eigenvalue( lyr[k_old].grid1d[i_old].qptr_h->ret_eval(E),
							   lyr[k_old].grid1d[i_old].qptr_h->ret_wave(E) )) {
		lyr[k].grid1d[i].qptr_h->modify_wave( lyr[k_old].grid1d[i_old].qptr_h->ret_eval(E),
						      lyr[k_old].grid1d[i_old].qptr_h->ret_wave(E) );
		if(fp) printf("WAVEFUNCTION_P:2 Eigenvalue already present\n");
	      }
	    }
	    if(lyr[k].grid1d[i].qptr_h->calc_prob(E,0)>lyr[k_old].grid1d[i_old].qptr_h->calc_prob(E,0)) d = secure_length ;
	    k_old = k;
	    i_old = i;
	    incr(k,i);
	    d += lyr[k].grid1d[i].xsize();
	  }
	}
      }
      
      if(!hole_wavefunction_normalization(E,fp)) nrerror("Wavefunction_p:: Not normalized!");
    }
  }



  // ----------------------------------------------------------------------
  if(switch4) {
    if(fp) printf("Find the lowest energy wavefunction\n");
    k = 0;
    i = 0;
    E = 1e10;
    while(test_upper_limit(k,i)) {
      if(lyr[k].grid1d[i].qptr_h!=0 && lyr[k].grid1d[i].qptr_h->ret_num()>0 && lyr[k].grid1d[i].qptr_h->ret_energy(0)<E )
	E = lyr[k].grid1d[i].qptr_h->ret_energy(0);
      incr(k,i);
    }

    
    if(E < 1e10) {
      if(fp) printf("Spot the number of peaks of the wave function\n");
      peak_count = 0;
      j = 0;
      k = 0;
      i = 0;
      prob1 = -1.0;
      prob2 = -1.0;
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_h==0 ) incr(k,i);
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_h!=0 && lyr[k].grid1d[i].qptr_h->check_norm(E)==-1 ) incr(k,i);
      if(lyr[k].grid1d[i].qptr_h!=0 && lyr[k].grid1d[i].qptr_h->check_norm(E)!=-1) prob1=lyr[k].grid1d[i].qptr_h->calc_prob(E,0);
      incr(k,i);
      if(lyr[k].grid1d[i].qptr_h!=0 && lyr[k].grid1d[i].qptr_h->check_norm(E)!=-1)  prob2=lyr[k].grid1d[i].qptr_h->calc_prob(E,0);
      incr(k,i);
      while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_h!=0 && lyr[k].grid1d[i].qptr_h->check_norm(E)!=-1 && prob1!=-1.0 && prob2!=-2.0) {
	prob3=lyr[k].grid1d[i].qptr_h->calc_prob(E,0);
	if(prob1<prob2 && prob2>prob3) peak_count++;
	prob1=prob2;
	prob2=prob3;
	incr(k,i);
      }
      if(fp) printf(" .... %d peaks found\n");
    }
    else {
      if(fp) printf("It seems there is not a single wavefunction.\n");
      peak_count = 2;      
    }

    if(fp) printf("Find the bottom of valence band\n");
    j = 1;
    j_pos = 0;
    E_min = 999.0;
    while(j<step_num-1) {
      if(band[j].energy<band[j-1].energy && band[j].energy<band[j+1].energy && band[j].energy<E_min) {
	E_min = band[j].energy;
	j_pos = j;
      }
      j++;
    }
      
    if(peak_count>1 && j_pos>0) {

      if(E==1e10) {
	j = 0;
	while(band[j_pos-j].energy<band[j_pos-j-1].energy && band[j_pos+j].energy<band[j_pos+j+1].energy) j++;
	E_max = (band[j_pos-j].energy>band[j_pos+j+1].energy)?band[j_pos+j].energy:band[j_pos-j].energy;
      }
      else E_max = E;
      if(E_min>=E_max) nrerror("WAVEFUNCTION_P: first wavefunction not found!");
      
      h = (E_max - E_min)*1e-4;
      E = E_min+h;
      found = 0;
      old_sum = 1e250;
      
      while(E<E_max && found !=1 ) {
	
	if(fp) printf("It creates and initilizes every Qntm_var class for each mesh point of the heterostructure\n"); 
	ptr_qntm = new Qntm_var * [xnum_mesh_points];
	tmp = new Eigenvalue;
	tmp->energy = E;
	tmp->gamma = 0;
	tmp->sigma = 0;
	tmp->ptr = 0;
	t = 0;
	k = 0;
	i = 0;
	while(test_upper_limit(k,i)) {
	  if(lyr[k].grid1d[i].qptr_h!=0) {
	    ptr_qntm[t] = new Qntm_var;
	    ptr_qntm[t]->set_eigen(tmp, lyr[k].grid1d[i].qptr_h->ret_Eo(),(lyr[k].grid1d[i].qptr_h->ret_offset()+lyr[k].grid1d[i].qptr_h->ret_leftedge())/LNTHCONV,lyr[k].grid1d[i].qptr_h->ret_leftedge()/LNTHCONV,lyr[k].grid1d[i].mh());
	  }
	  else ptr_qntm[t] = 0;
	  incr(k,i);
	  t++;
	}	  

	if(fp) printf("Set the complex coefficients of the wave function\n");
	bb = 1.0;
	sum = 0.0;
	k = 0;
	i = 0;
	t = 0;
	while(test_upper_limit(k,i) && ptr_qntm[t]==0 ) {
	  incr(k,i);
	  t++;
	}
	if(test_upper_limit(k,i)) {
	  Ac.rl = 0.0;
	  Ac.im = 0.0;
	  Bc.rl = bb;
	  Bc.im = 0.0;
	  if(!ptr_qntm[t]->set_wave(E,Ac,Bc)) nrerror("H_NORMALIZTION out of sync");
	  sum += ptr_qntm[t]->calc_prob(E,lyr[k].grid1d[i].xsize());
	  incr(k,i);
	  t++;
	  while( ptr_qntm[t]!=0 && ptr_qntm[t]->check_norm(E)!=-1) {
	    if(fabs(ptr_qntm[t]->ret_Eo()-ptr_qntm[t-1]->ret_Eo())>1e-6)
	      ptr_qntm[t]->propagation(E,Ac,Bc,ptr_qntm[t-1]);
	    ptr_qntm[t]->set_wave(E,Ac,Bc);
	    sum += ptr_qntm[t]->calc_prob(E,lyr[k].grid1d[i].xsize());
	    incr(k,i);
	    t++;
	  }
	}
	
	if(fp) printf("Spotting the number of peaks\n");
	n_max = 0;
	t = 1;
	while(t<xnum_mesh_points-1) {
	  ckr = 0;
	  if(ptr_qntm[t-1] && ptr_qntm[t] && ptr_qntm[t+1]) {
	    if(ptr_qntm[t]->calc_prob(E,0)>ptr_qntm[t-1]->calc_prob(E,0) 
	       && ptr_qntm[t]->calc_prob(E,0)>ptr_qntm[t+1]->calc_prob(E,0)) ckr = 1;
	  }
	  if(ckr) n_max++;
	  t++;
	}
	if(fp) printf("In E=%le I found %d peaks\n",E,n_max);

	if(n_max==1 && sum>old_sum) {
	  fopen_s(&fp1,"OUTPUT_FILES/tempo_p.dat","w");
	  t = 0;
	  while(t<xnum_mesh_points) {
	    if(ptr_qntm[t]) fprintf(fp1,"%d  %le\n",t,ptr_qntm[t]->calc_prob(E,0));
	    t++;
	  }
	  fclose(fp1);
	  found = 1;
	  E -= h;
	}
	else {
	  old_sum = sum;
	  E += h;
	}

	if(fp) printf("Removing temporary arrays\n"); 
	delete tmp;
	for(t=0;t<xnum_mesh_points;t++) {
	  if(ptr_qntm[t]!=0) delete ptr_qntm[t];
	  ptr_qntm[t] = 0;
	}
	delete [] ptr_qntm;
	ptr_qntm = 0;
	
      }
      
      if(found) {
	
	j = 0;
	k = 0;
	i = 0;
	t = 0;
	while(test_upper_limit(k,i)) {
	  if(lyr[k].grid1d[i].qptr_h) lyr[k].grid1d[i].qptr_h->add_eigenvalue(E,lyr[k].grid1d[i].mh());
	  incr(k,i);
	}
	if(fp) printf("Normalization of the wave function\n");
	if(!hole_wavefunction_normalization(E,fp)) nrerror("Wavefunction_p:: Not normalized!");
	
	for(ll=ndone-1;ll>=0;ll--) done[ll+1]=done[ll];
	  done[0]=E;
	  ndone++;
      }
      else printf("WAVEFUNCTION_P: Lowest value not found!\n");
      
    }
  }


  if(fp) printf("Set the probability vector and dump it in a file\n");
  for(ll=0;ll<ndone;ll++) {    
    E = done[ll];
    fopen_s(&fp1,str1,"w");
    j = 0;
    k = band[0].left[0];
    i = band[0].left[1];
    while(j<step_num) {
      while(test_upper_limit(k,i,band[j].right[0],band[j].right[1]) ) {
	lyr[k].grid1d[i].qptr_h->set_prob_h(E,lyr[k].grid1d[i].xsize(),fp1);
	incr(k,i);
      }
      j++;
    }
    fclose(fp1);
    str1[19]++;
  }


}



/* ------------------------------------------------ */
/* ------------------------------------------------ */


Eigenvalue ** Device::eigenvalues(int Number_of_Steps, Step * Step_Band_Profile, FILE * fp, FILE * fp1) {

  short switch1 = 1;
  int i,k,num_pks,num_pks1,num_En,tail;
  const int MAX_PK = 50;
  double E1,E2,** tras,E_STEP = 1e-3;
  Eigenvalue tmp_eigen, * tmp_pnt_eigen,* eval0,* eval1, ** store;
  FILE * fp2;

  int Number_of_eigenvalues_required = sim->Number_of_eigenvalues;


  if(fp) printf("EIGENVALUES2 STARTING ....");

  if( Number_of_eigenvalues_required>MAX_PK) {
    printf("You required more eigenvalue than the program is able to cope with\n");
    printf("The limit is set at %d eigenvalue\n",MAX_PK);
    printf("If you want to increase the limit you must modify the MAX_PK constant in file qdevice.cpp function eigenvalues2() and recompile.");
  }

  if(fp) printf("Allocation of temporary 2D array\n");
  store = new Eigenvalue* [Number_of_Steps];
  if(store == NULL) nrerror("Dynamic allocation failure of store in eigenvalues");
  for(i=0;i<Number_of_Steps;i++) store[i]=0;


  if(switch1) {
    if(fp) printf("Call to left_to_right Transmission coefficient function.\n");
    tras = transmission_coefficient(Number_of_Steps,Step_Band_Profile,0,0,-2,num_En,fp,0,0);
    if(fp) {
      fopen_s(&fp2,"OUTPUT_FILES/TT.dat","w");
      for(i=0;i<num_En;i++) fprintf(fp2,"%lf  %le\n",tras[0][i],tras[1][i]);
      fclose(fp2);
    }
  }
  else {
    if(fp) printf("Call to Right_to_left Transmission coefficient function.\n");
    tras = transmission_coefficient_bidir(Number_of_Steps,Step_Band_Profile,0,0,-2,num_En,fp,-1,-1);
    if(fp) {
      fopen_s(&fp2,"OUTPUT_FILES/TTR.dat","w");
      for(i=0;i<num_En;i++) fprintf(fp2,"%lf  %le\n",tras[0][i],tras[1][i]);
      fclose(fp2);
    }    
  }

  if(fp) printf("Call to Lorenzian fitting function. It returns number and position of the peaks in the trasmission coefficient.\n");
  tail = 10; // It gives the resolution used for spotting a peak in the lorenzian fitting function. The higher, the less noise is spotted! 
  eval0 = lorenzian( num_En, tras[1], tras[0], num_pks, tail  , Number_of_eigenvalues_required );
  
  delete [] tras[0];
  delete [] tras[1];
  delete [] tras;
  tras = 0;

  if(fp) printf("FOR EACH PEAK WE REPEAT ITS CHARACTERIZATION MORE PRECISELY\n");
  k=0;
  while(k<num_pks) {

    E1 = eval0[k].energy - E_STEP;
    E2 = eval0[k].energy + E_STEP;

    tras = transmission_coefficient(Number_of_Steps,Step_Band_Profile,E1,E2,2.0*E_STEP*E_STEP,num_En,fp,0,0);
    if(fp) {
      fopen_s(&fp2,"OUTPUT_FILES/TT0.dat","w");
      for(i=0;i<num_En;i++) fprintf(fp2,"%lf  %lf\n",tras[0][i],tras[1][i]);
      fclose(fp2);
    }

    tail = 100;
    eval1 = lorenzian( num_En, tras[1], tras[0], num_pks1, tail, 1);

    while( num_pks1==0 && tail>10 ) {
      tail -= 2;
      delete [] eval1;
      eval1 = 0;
      eval1 = lorenzian( num_En,tras[1],tras[0],num_pks1,tail,1);
    }

    if(num_pks1==0) {
      if(fp) {
	printf("CALC_EIGENVALUES: %d peaks obtained, 1 expected!\n",num_pks1);
	printf("num=%d at Energy=%le\n",k,eval0[k].energy);
	fopen_s(&fp2,"OUTPUT_FILES/err_eigen.dat","w");
	for(i=0;i<num_En;i++) fprintf(fp2,"%lf  %lf  \n",tras[0][i],tras[1][i]);
	fclose(fp2);
      }
    }
    else eval0[k] = eval1[0];
    k++;

    delete [] tras[0];
    delete [] tras[1];
    delete [] tras;
    tras = 0;
    delete [] eval1;
    eval1 = 0;
  }


  if(fp) printf("Allocate the 2D array containing for each step of the band structure a list of eigenvalues");
  if(num_pks>0) {
    for(i=0;i<Number_of_Steps;i++) {
      if(store[i]!=0) {
	tmp_pnt_eigen = store[i];
	while(tmp_pnt_eigen->ptr!=NULL) tmp_pnt_eigen = tmp_pnt_eigen->ptr;
	for(k=0;k<num_pks;k++) {
	  tmp_pnt_eigen->ptr = new Eigenvalue;
	  *(tmp_pnt_eigen->ptr) = eval0[k];
	  tmp_pnt_eigen->ptr->ptr = 0;
	  tmp_pnt_eigen = tmp_pnt_eigen->ptr;
	}
      }
      else {
	store[i] = new Eigenvalue;
	*store[i] = eval0[0];
	tmp_pnt_eigen = store[i];
	for(k=1;k<num_pks;k++) {
	  tmp_pnt_eigen->ptr = new Eigenvalue;
	  *(tmp_pnt_eigen->ptr) = eval0[k];
	  tmp_pnt_eigen->ptr->ptr = 0;
	  tmp_pnt_eigen = tmp_pnt_eigen->ptr;
	}
      }
    }
  }
  delete [] eval0;

  if(fp) printf("EIGENVALUES2 ENDS");

  return store;


}

/* -------------------------------------------------- */
/* -------------------------------------------------- */

double ** Device::transmission_coefficient_bidir(int num_step,Step * band,double E1,double E2,double gap,int & num_en,FILE * fp,int initial,int final)	{

  int i,j,num,dir;
  const int MIN = 500;
  const double OVER_TOP = 1.0; 
  double d,E_max,E_min,E,E_fin,E_in,E_gap;

  double ** tr = 0;
  Complex * k_vec = 0,prod;
  Tmatrix * Tmat = 0;


  if(final>initial) {
    num = final - initial +1;
    dir = 1;
  }
  else if(initial>final) {
    num = initial - final +1;
    dir = -1;
  }
  else if(initial==final && initial>=0) {
    initial = 0;
    num = num_step;
    dir = 1;
  }
  else if(initial==final && initial<0){
    initial = num_step-1;
    num = num_step;
    dir = -1;
  }

  if(fp) fprintf(fp,"Calculation of the transmission coefficient ...\n");

  E_gap = gap;
  E_in 	= E1;
  E_fin = E2;

  k_vec	= new Complex [num];
  if(k_vec==NULL) nrerror("Memory allocation failure for k_vec in transmission_coefficient");
  Tmat = new Tmatrix [num-1];
  if(Tmat==NULL) nrerror("Memory allocation failure for Tmat in transmission_coefficient");

  /* SET THE TYPE OF CALCULATION */
  if(E_gap<0) {
    E_max = -50.0;
    E_min = 50.0;
    for(j=0;j<num;j++) {
      if(E_max<band[j*dir+initial].energy)	E_max=band[j*dir+initial].energy;
      if(E_min>band[j*dir+initial].energy)	E_min=band[j*dir+initial].energy;
    }
    E_in  = E_min;
    E_fin = E_max;
    num_en= MIN;
    E_gap = (E_fin - E_in )/(double)MIN;
  }
  else if(E_gap==0) {
    E_max = -50.0;
    for(j=0;j<num;j++) if(E_max<band[j*dir+initial].energy) E_max=band[j*dir+initial].energy;
    E_in  = E_max;
    E_fin = E_max+OVER_TOP;
    num_en= MIN;
    E_gap = OVER_TOP/(double)MIN;
  }
  else if(E_fin==E_in) num_en = 1;
  else num_en = int((double)(E_fin-E_in )/E_gap);

  /* Calculation at least over MIN points */
  if(num_en<MIN) {
    num_en = MIN;
    E_gap  = (E_fin - E_in )/(double)MIN;
  }

  /* ALLOCATION ARRAYS */
  tr = new double* [2];
  if(tr==NULL) nrerror("Memory allocation failure for tr in transmission_coefficient");
  tr[0] = new double [num_en];
  if(tr[0]==NULL) nrerror("Memory allocation failure for tr[0] in transmission_coefficient");
  tr[1] = new double [num_en];
  if(tr[1]==NULL) nrerror("Memory allocation failure for tr[1] in transmission_coefficient");

  E = E_in + 0.1*E_gap;
  for(i=0;i<num_en;i++)	{ // Iteration on Energy
    for(j=0;j<num;j++) { // Iteration on step.
      if(fabs(E-band[j*dir+initial].energy)<1.0e-20) {
	printf("Transmission_coefficient: Energy equals to step: correction made!\n");
	E = E+0.01*E_gap;
      }
      k_vec[j]=cplx_sqrt(band[j*dir+initial].mass*KKK*(E-band[j*dir+initial].energy));
    }
    d =	0.0;
    for(j=0;j<num-1;j++)	{
      Tmat[j].is(k_vec[j],k_vec[j+1],band[j*dir+initial].mass,band[(j+1)*dir+initial].mass,d);
      d	+= band[(j+1)*dir+initial].size*LNTHCONV;
    }
			
    prod = product(num-1,Tmat);

    if(prod.rl==0 && prod.im==0) {
      /* printf("CALC_TRANSMISSION: Zero transmission coefficient at %d  E=%le in [%le,%le]\n",i,E, E_in,E_fin);
	 delete [] tr[0];
	 delete [] tr[1];
	 delete [] tr;
	 delete [] k_vec;
	 delete [] Tmat;
	 exit(1);
      */
      if(i>0)  tr[1][i]=tr[1][i-1];
      else tr[1][i]=1e-20;
    }
    else tr[1][i]=cplx_mod_sqr(prod);

    tr[0][i]=E;
    E += E_gap;

  }

  delete [] k_vec;
  delete [] Tmat;

  if(fp) fprintf(fp,"Transmission coefficient done!\n");

  return tr;

}

/* ----------------------------- */
/* ----------------------------- */

double ** Device::optical_transitions(int & num_transition) {

  
  int i,k,j;
  double transition[MAX_TRANSITIONS],overlap[MAX_TRANSITIONS];
  double tmp,** result,E,d,d0,h,x,E_e,E_h,off1,off2;
  Fwave w1,w2;

  /* SPOT THE TRANSITIONS */
  num_transition = 0;
  k = 0;
  i = 0;
  while(k<num_lyr-1 || (k==num_lyr-1 &&  i<lyr[k].numx())) {
    if(lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_h!=0) {
      if(lyr[k].grid1d[i].qptr_e->ret_num()>0 &&  lyr[k].grid1d[i].qptr_h->ret_num()>0) {
	tmp = lyr[k].grid1d[i].qptr_e->ret_energy(0)+lyr[k].grid1d[i].qptr_h->ret_energy(0);
	if( num_transition==0 || tmp!=transition[num_transition-1] ) {
	  transition[num_transition] = tmp;
	  num_transition++;
	}
      }
    }
    if(i<lyr[k].numx()-1) i++;
    else { k++; i=0; }
  }

  /* CALCULATE THE OVERLAP */
  for(j=0;j<num_transition;j++) {
    k = 0;
    i = 0;
    d = 0;
    overlap[j] = 0;
    E = transition[j];
    while(k<num_lyr-1 || (k==num_lyr-1 &&  i<lyr[k].numx())) {
      d0 = d;
      d += lyr[k].grid1d[i].xsize();
      if(lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_h!=0) {
	if(lyr[k].grid1d[i].qptr_e->ret_num()>0 &&  lyr[k].grid1d[i].qptr_h->ret_num()>0) {
	  E_e = lyr[k].grid1d[i].qptr_e->ret_energy(0);
	  E_h = lyr[k].grid1d[i].qptr_h->ret_energy(0);
	  if(fabs((E_e+E_h)-E)<1e-5) {
	    x = d0*LNTHCONV;
	    h = LNTHCONV*(d - d0)/50.0;
	    while(x<(d*LNTHCONV)) {
	      w1 = lyr[k].grid1d[i].qptr_e->ret_wave(E_e);
	      w2 = lyr[k].grid1d[i].qptr_h->ret_wave(E_h);
	      off1 = lyr[k].grid1d[i].qptr_e->ret_leftedge();
	      off2 = lyr[k].grid1d[i].qptr_h->ret_leftedge();
	      overlap[j] += h*cplx_mod_sqr(w1.eval(x-off1)*w2.eval(x-off2));
	      x += h;
	    }
	  }
	  else {
	    printf("E=%lf  E_e+E_h=%lf\n",E,E_e+E_h);
	    nrerror("OPTICAL_TRANSITION: No matching energies");
	  }
	}
      }
      if(i<lyr[k].numx()-1) i++;
      else { k++; i=0; }
    }
  }

  result = new double * [2];
  result[0] = new double [num_transition];
  result[1] = new double [num_transition];
  for(i=0;i<num_transition;i++) {
    result[0][i]=transition[i];
    result[1][i]=overlap[i];
  }
  return result;

}




/* ----------------------------- */
/* ----------------------------- */

double Device::field_in_qw() {

  int k,i;

  k=0;
  i=0;
  while(test_upper_limit(k,i) && (lyr[k].grid1d[i].qptr_e==0 || lyr[k].grid1d[i].qptr_h==0) ) incr(k,i);
  if(!test_upper_limit(k,i)) return 0;
  while(test_upper_limit(k,i) && i!=10) incr(k,i);
  return lyr[k].grid1d[i].electric_field();

}

/* ----------------------------- */
/* ----------------------------- */

int Device::hole_wavefunction_normalization(double E, FILE * fp) {

  int i,k,i_old,k_old,count;
  double bb,sum;
  Complex Ac,Bc;

  bb = 1.0;
  do {
    count = 0;
    sum	= 0.0;
    k = 0;
    i = 0;
    while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_h==0 ) incr(k,i);
    if(!test_upper_limit(k,i)) return 0;
    Ac.rl = 0.0;
    Ac.im = 0.0;
    Bc.rl = bb;
    Bc.im = 0.0;
    if(!lyr[k].grid1d[i].qptr_h->set_wave(E,Ac,Bc)) {
      nrerror("H_NORMALIZTION out of sync");
      sum += lyr[k].grid1d[i].qptr_h->calc_prob(E,lyr[k].grid1d[i].xsize());
    }
    k_old = k;
    i_old = i;
    incr(k,i);
    while( lyr[k].grid1d[i].qptr_h!=0 && lyr[k].grid1d[i].qptr_h->check_norm(E)!=-1) {
      if(fabs(lyr[k].grid1d[i].qptr_h->ret_Eo()-lyr[k_old].grid1d[i_old].qptr_h->ret_Eo())>1e-6)
	lyr[k].grid1d[i].qptr_h->propagation(E,Ac,Bc,lyr[k_old].grid1d[i_old].qptr_h);
      lyr[k].grid1d[i].qptr_h->set_wave(E,Ac,Bc);
      sum += lyr[k].grid1d[i].qptr_h->calc_prob(E,lyr[k].grid1d[i].xsize());
      k_old = k;
      i_old = i;
      incr(k,i);
    }
    count++;
    if(sum!=0) bb = bb/sqrt(sum);
    else nrerror("Ohps in hole_wavefunction_normalization: sum ==0\n");
    if(fp) printf("count=%d, sum=%le\n",count,sum);
  } while(fabs(sum-1.0)>1e-3 && count<1000);

  if(count<1000) return 1;
  else return 0;
}

/* ----------------------------- */
/* ----------------------------- */

int Device::electron_wavefunction_normalization(double E, FILE * fp) {

  int i,k,i_old,k_old,count;
  double bb,sum;
  Complex Ac,Bc;

  bb = 1.0;
  do {
    count = 0;
    sum	= 0.0;
    k = 0;
    i = 0;
    while(test_upper_limit(k,i) && lyr[k].grid1d[i].qptr_e==0 ) incr(k,i);
    if(!test_upper_limit(k,i)) return 0;
    Ac.rl = 0.0;
    Ac.im = 0.0;
    Bc.rl = bb;
    Bc.im = 0.0;
    if(!lyr[k].grid1d[i].qptr_e->set_wave(E,Ac,Bc)) {
      nrerror("N_NORMALIZTION out of sync");
      sum += lyr[k].grid1d[i].qptr_e->calc_prob(E,lyr[k].grid1d[i].xsize());
    }
    k_old = k;
    i_old = i;
    incr(k,i);
    while( lyr[k].grid1d[i].qptr_e!=0 && lyr[k].grid1d[i].qptr_e->check_norm(E)!=-1) {
      if(fabs(lyr[k].grid1d[i].qptr_e->ret_Eo()-lyr[k_old].grid1d[i_old].qptr_e->ret_Eo())>1e-6)
	lyr[k].grid1d[i].qptr_e->propagation(E,Ac,Bc,lyr[k_old].grid1d[i_old].qptr_e);
      lyr[k].grid1d[i].qptr_e->set_wave(E,Ac,Bc);
      sum += lyr[k].grid1d[i].qptr_e->calc_prob(E,lyr[k].grid1d[i].xsize());
      k_old = k;
      i_old = i;
      incr(k,i);
    }
    count++;
    if(sum!=0) bb = bb/sqrt(sum);
    else nrerror("Ohps in electron_wavefuntion_normalization: sum ==0\n");
    if(fp) printf("count=%d, sum=%le\n",count,sum);
  } while(fabs(sum-1.0)>1e-3 && count<1000);

  if(count<1000) return 1;
  else return 0;
}



/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */

void Device::Schrodinger(short debug) {

  int j,i,**left=0,**right=0,num_min,num_step_cb,num_step_vb,num_en;
  double ** trans=0;
  Step * cb_stp=0;
  Step * vb_stp=0;
  Eigenvalue ** eeval =0,** heval =0, * eig_pnt =0;
  FILE * fp, *fp1;

  if(debug) printf("SCHRODINGER SOLVER ...\n");

  char namefile[30]="OUTPUT_FILES/e_tras0.dat";
  char namefile2[30]="OUTPUT_FILES/cb_step0.dat";
  char namefile3[30]="OUTPUT_FILES/cb_eval0.dat";
  char namefile4[30]="OUTPUT_FILES/nwave00.dat";
  char namefile5[30]="OUTPUT_FILES/pwave00.dat";
  char file1[30]="OUTPUT_FILES/n_even.dat";
  char file2[30]="OUTPUT_FILES/n_odd.dat";
  char file3[30]="OUTPUT_FILES/n_even.dat";
  char file4[30]="OUTPUT_FILES/n_odd.dat";

  if(debug) fopen_s(&fp,"OUTPUT_FILES/quantum.log","w");
  else fp = 0;

  q_eraser(); /* Destroy the previous quantum structure */


  // SOLVE ANALYTICALLY 
  //printf("V=%le  x=%le  me0=%le  meQW=%le\n",fabs(lyr[1].grid1d[lyr[0].numx()-1].cb()-lyr[2].grid1d[0].cb()),lyr[2].xsize()*LNTHCONV,lyr[1].grid1d[lyr[0].numx()-1].me(),lyr[2].grid1d[0].me());
  //analytical(fabs(lyr[1].grid1d[lyr[0].numx()-1].cb()-lyr[2].grid1d[0].cb()),10.0,lyr[2].xsize()*LNTHCONV,lyr[1].grid1d[lyr[0].numx()-1].me(),lyr[2].grid1d[0].me(),file1,file2);


  // CONDUCTION BAND
  left 	= new int* [2];
  right = new int* [2];
  if(debug) printf("--->find_hetero_cb...\n");
  num_min = find_hetero_cb(left,right,fp);
  if(debug) printf("--->...(%d minima)!\n",num_min);


  for(j=0;j<num_min;j++) { /* Iteration on quantum structures */

    if(debug) printf("------>cb=%d : (%d,%d)<-->(%d,%d)\n",j,left[0][j],left[1][j],right[0][j],right[1][j]);

    /* STEP APPROXIMATION */
    if(debug) printf("------>Step_approximation_cb...");
    fopen_s(&fp1,namefile2,"w");
    cb_stp = step_approximation_cb(num_step_cb,left[0][j],left[1][j],right[0][j],right[1][j],fp,fp1);
    fclose(fp1);
    if(debug) printf("... %d num_step in CB\n",num_step_cb);

    /* EIGENVALUE CALCULATION */
    if(debug) printf("------>Eigenvalues...");
    fopen_s(&fp1,namefile3,"w");
    eeval = eigenvalues(num_step_cb,cb_stp,fp,fp1);
    fclose(fp1);
    
    /* CALCULATE THE WAVEFUNCTIONS */
    if(debug) printf("------>Wavefunction_n...\n");
    namefile4[19]='0';
    wavefunction_n(eeval,num_step_cb,cb_stp,fp,namefile4);
    delete [] eeval;

    /* CALCULATE TRANSMISSION COEFFICIENT */
    if(debug) {
      printf("------>Transmission_coefficient...\n");
      trans = transmission_coefficient(num_step_cb,cb_stp,0,0,0,num_en,fp);
      fopen_s(&fp1,namefile,"w");
      for(i=0;i<num_en;i++) fprintf(fp1,"%le  %le\n",trans[0][i],trans[1][i]);
      fclose(fp1);
      if(debug) printf("------>File %s dumped\n",namefile);
      namefile[19]++;
      if(trans) {
	delete [] trans[0];
	delete [] trans[1];
	delete [] trans;
      }
    }

    namefile2[20]++;
    namefile3[20]++;
    namefile4[18]++;

    if(cb_stp) delete [] cb_stp;
    if(left) {
      delete [] left[0];
      delete [] left[1];
      delete [] left;
    }
    if(right) {
      delete [] right[0];
      delete [] right[1];
      delete [] right;
    }
  }
  
  /* VALENCE BAND */
    
  namefile[13]='h';
  namefile[19]='0';
  namefile2[13]='v';
  namefile2[20]='0';
  namefile3[13]='v';
  namefile3[20]='0';

  left 	= new int* [2];
  right = new int* [2];
  if(debug) printf("--->find_hetero_vb...\n");
  num_min = find_hetero_vb(left, right,fp);
  if(debug) printf("--->...(%d minima)!\n",num_min);

  for(j=0;j<num_min;j++) {

    if(debug) printf("------>vb=%d : (%d,%d)<-->(%d,%d)\n",j,left[0][j],left[1][j],right[0][j],right[1][j]);

    if(debug) printf("------>Step_approximation_vb...\n");
    fopen_s(&fp1,namefile2,"w");
    vb_stp = step_approximation_vb(num_step_vb,left[0][j],left[1][j],right[0][j],right[1][j],fp,fp1);
    fclose(fp1);

    if(debug) printf("------>Eigenvalues...\n");
    fopen_s(&fp1,namefile3,"w");
    heval = eigenvalues(num_step_vb,vb_stp,fp,fp1);
    fclose(fp1);

    if(debug) printf("--->Wavefunction_p...\n");
    namefile5[19]='0';
    wavefunction_p(heval,num_step_vb,vb_stp,fp,namefile5);
    delete [] heval;

    if(debug) {
      printf("--->VB Transmission_coefficient...\n");
      trans = transmission_coefficient( num_step_vb,vb_stp,0,0,0,num_en,fp);
      fopen_s(&fp1,namefile,"w");
      for(i=0;i<num_en;i++) fprintf(fp1,"%le  %le\n",trans[0][i],trans[1][i]);
      fclose(fp1);
      namefile[19]++;
      if(trans) {
	delete [] trans[0];
	delete [] trans[1];
	delete [] trans;
      }
    }

    namefile2[20]++;
    namefile5[18]++;

    if(vb_stp) delete [] vb_stp;
    if(left) {
      delete [] left[0];
      delete [] left[1];
      delete [] left;
    }
    if(right) {
      delete [] right[0];
      delete [] right[1];
      delete [] right;
    }
    
  }



  if(debug) fclose(fp);
  if(debug) printf("...DONE!\n\n");
}
